/*
 * modelcodon.cpp
 *
 *  Created on: May 24, 2013
 *      Author: minh
 */

#include "modelcodon.h"
#include <string>


const double MIN_OMEGA_KAPPA = 0.001;
const double MAX_OMEGA_KAPPA = 50.0;

/* Empirical codon model restricted (Kosiol et al. 2007), source: http://www.ebi.ac.uk/goldman/ECM/ */
string model_ECMrest1 =
"11.192024 \
1.315610 0.010896 \
5.427076 4.756288 24.748755 \
1.658051 0.000000 0.000000 0.000000 \
0.000000 1.913571 0.000000 0.000000 13.889102 \
0.000000 0.000000 2.952332 0.000000 44.407955 13.681751 \
0.000000 0.000000 0.000000 8.126914 17.057443 65.097021 12.991861 \
6.610894 0.000000 0.000000 0.000000 2.206054 0.000000 0.000000 0.000000 \
0.000000 5.177930 0.000000 0.000000 0.000000 5.615472 0.000000 0.000000 19.942818 \
3.347364 0.000000 0.000000 0.000000 6.191481 0.000000 0.000000 0.000000 0.582084 0.000000 \
0.000000 1.558523 0.000000 0.000000 0.000000 9.339206 0.000000 0.000000 0.000000 0.144278 44.777964 \
0.000000 0.000000 0.000000 5.369644 0.000000 0.000000 0.000000 4.662001 0.000000 0.000000 0.677177 0.073268 \
2.090751 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
0.000000 2.266373 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.905484 \
0.000000 0.000000 75.752638 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 56.803876 7.811205 \
0.000000 0.000000 0.000000 20.877218 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.432339 22.078564 5.650116 \
0.000000 0.000000 0.000000 0.000000 1.769355 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.263838 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 2.704601 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.389735 0.000000 0.000000 17.461627 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.312811 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.393680 0.000000 35.480963 12.053827 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.303480 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.477616 8.407091 28.557939 11.295213 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.444964 0.000000 0.000000 0.000000 0.000000 1.583116 0.000000 0.000000 0.000000 1.021682 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.087801 0.000000 0.000000 0.000000 0.000000 3.230751 0.000000 0.000000 0.000000 3.774544 0.000000 0.000000 28.086160 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.419058 0.000000 0.000000 0.000000 5.381868 0.000000 3.440380 1.918904 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.812540 0.000000 0.000000 0.000000 1.794388 1.086327 5.369463 14.959151 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.617091 0.000000 0.000000 0.779565 0.000000 0.000000 0.000000 0.334165 0.000000 0.000000 0.000000 3.019726 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.632945 0.000000 0.000000 2.250770 0.000000 0.000000 0.000000 1.699302 0.000000 0.000000 0.000000 7.016899 0.000000 0.000000 14.603857 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.023939 0.000000 0.000000 0.000000 1.693662 0.000000 0.000000 0.000000 6.415757 0.000000 99.459951 14.930266 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.026086 0.000000 0.000000 0.000000 1.462945 0.000000 0.000000 0.000000 3.144296 0.000000 0.000000 0.000000 19.920977 30.804750 79.483730 13.919752 \
1.682029 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.301225 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
0.000000 0.786043 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.381841 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.140728 \
0.000000 0.000000 10.116588 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.134459 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 18.298900 4.623936 \
0.000000 0.000000 0.000000 7.911096 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.570123 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.281784 1.303951 2.082128 \
0.000000 0.000000 0.000000 0.000000 38.229100 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.578976 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.801564 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 15.793595 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.434550 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.231468 0.000000 0.000000 6.035740 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.033932 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.925575 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.962350 0.000000 28.307876 6.967655 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 17.103904 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.238450 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.155285 19.578982 38.414969 12.678802 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.245405 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.004762 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.501054 0.000000 0.000000 0.000000 11.715476 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.228361 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.105602 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.292691 0.000000 0.000000 0.000000 2.134740 0.000000 0.000000 13.863648 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.404436 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.647620 0.000000 0.000000 0.000000 3.919360 0.000000 4.929483 0.366267 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.715692 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.975074 0.000000 0.000000 0.000000 5.869857 1.010212 0.982893 10.762877 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.719489 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.834666 0.000000 0.000000 0.000000 0.578118 0.000000 0.000000 0.000000 39.399322 0.000000 0.000000 0.000000 16.623529 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.047654 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.033630 0.000000 0.000000 0.000000 0.437779 0.000000 0.000000 0.000000 21.337943 0.000000 0.000000 0.000000 7.784768 0.000000 0.000000 26.637668 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 92.372238 0.000000 0.000000 0.000000 1.903175 0.000000 0.000000 0.000000 0.754055 0.000000 0.000000 0.000000 8.423762 0.000000 1.792245 0.120900 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.825082 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 133.296291 0.000000 0.000000 0.000000 2.231662 0.000000 0.000000 0.000000 22.577271 0.000000 0.000000 0.000000 21.000358 3.324581 6.011970 36.292705 \
2.261813 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.473623 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.096281 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
0.000000 1.923392 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.914972 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.137337 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.669955 \
0.000000 0.000000 2.362720 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.737489 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 25.294298 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 26.045078 3.531461 \
0.000000 0.000000 0.000000 2.022101 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.164805 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.078444 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.901167 21.657664 11.898141 \
0.000000 0.000000 0.000000 0.000000 5.540052 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.159185 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.107629 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.682092 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 7.675838 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.120189 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.312255 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.308415 0.000000 0.000000 6.516319 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 9.880382 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.923972 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.064069 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.291148 0.000000 21.910225 5.090423 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 21.863158 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.034856 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 25.461549 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.166554 5.512586 20.715347 9.529141 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.367553 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.383706 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.091654 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.352915 0.000000 0.000000 0.000000 0.693026 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.294702 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.006827 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.686074 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.208522 0.000000 0.000000 0.000000 1.866565 0.000000 0.000000 10.605899 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.485369 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.811398 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.277861 0.000000 0.000000 0.000000 2.774445 0.000000 2.710610 0.650088 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.686782 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.090641 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.476105 0.000000 0.000000 0.000000 9.441919 1.296294 3.779053 10.153570 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.104727 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.041150 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.590780 0.000000 0.000000 0.000000 0.503385 0.000000 0.000000 0.000000 1.541379 0.000000 0.000000 0.000000 1.042624 0.000000 0.000000 0.000000 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.552851 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.252470 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.285543 0.000000 0.000000 0.000000 0.542717 0.000000 0.000000 0.000000 2.303487 0.000000 0.000000 0.000000 1.561629 0.000000 0.000000 9.488520";
string model_ECMrest = model_ECMrest1 + 
"0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.091041 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.432410 0.000000 0.000000 0.000000 0.702411 0.000000 0.000000 0.000000 2.985093 0.000000 0.000000 0.000000 0.874000 0.000000 20.518100 4.120953 \
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.810856 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.803738 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.388514 0.000000 0.000000 0.000000 0.302501 0.000000 0.000000 0.000000 6.644971 0.000000 0.000000 0.000000 1.393810 13.246936 18.064826 19.084271 \
\
0.022103  0.021383  0.016387  0.015425  0.011880  0.011131  0.009750  0.008956  0.015965  0.015782  0.006025  0.007029  0.011880  0.014467  0.017386  0.007600  0.028839  0.010007  0.010100  0.010642  0.011843  0.011097  0.011703  0.016076  0.020211  0.008311  0.014148  0.004800  0.007837  0.025576  0.023441  0.013551  0.020102  0.013424  0.020201  0.015528  0.012142  0.023006  0.020171  0.030001  0.026344  0.010142  0.011679  0.010372  0.008195  0.019047  0.018938  0.010901  0.022747  0.019005  0.028307  0.015908  0.018853  0.028198  0.024532  0.033223  0.031878  0.016852  0.022982  0.015796  0.010191 \
\
TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
GGG";

/* Empirical codon model unrestricted (Kosiol et al. 2007), source: http://www.ebi.ac.uk/goldman/ECM/ */
string model_ECMunrest1 =
"16.011531 \
2.395822 0.151858 \
1.204356 0.675537 18.541946 \
0.773935 0.052602 0.249707 0.274990 \
0.030074 0.656004 0.011609 0.158873 23.655090 \
0.278090 0.056677 1.184813 0.611887 35.921779 15.982573 \
0.034137 0.198277 0.010188 0.694091 11.510965 35.359077 17.424222 \
4.317981 0.503397 0.798582 0.337279 0.688169 0.047115 0.341791 0.058136 \
0.481042 4.483501 0.033529 0.177833 0.069588 0.524116 0.070809 0.213967 24.177765 \
0.733587 0.076912 0.645571 0.395942 1.811753 0.343463 0.751980 0.143447 0.822999 0.054860 \
0.045951 0.561620 0.040012 0.240632 0.138244 1.323765 0.121937 0.493179 0.068342 0.628438 56.838378 \
0.786871 1.183337 0.271072 0.632947 0.069758 0.081312 0.195833 0.410046 1.140051 1.421996 0.264556 0.210115 \
2.016257 0.207692 12.035723 11.161511 0.277929 0.000186 0.000289 0.000000 0.485469 0.000299 0.543240 0.000674 0.010122 \
0.083684 2.306110 1.373823 5.651603 0.000085 0.342813 0.000096 0.000344 0.000116 0.622089 0.000466 0.674176 0.113701 15.874441 \
1.036474 0.198558 27.219895 16.560966 0.000678 0.000186 0.496046 0.000115 0.016650 0.011978 0.020649 0.021578 0.017106 21.437257 8.808275 \
0.073550 1.341144 1.045943 12.455337 0.000000 0.001022 0.000000 0.266943 0.004815 0.308859 0.002639 0.265948 0.504866 4.802017 15.484088 8.319767 \
0.324368 0.000141 0.001358 0.003499 2.846677 0.196358 0.544474 0.078776 0.337879 0.000479 0.239715 0.000270 0.061833 0.822643 0.036254 0.181411 0.014388 \
0.000140 0.285635 0.000000 0.000382 0.101204 2.487136 0.072352 0.432520 0.000116 0.310416 0.000000 0.215779 0.032564 0.026571 0.648769 0.040087 0.149771 23.496083 \
0.025217 0.006558 0.261069 0.005535 0.487542 0.138742 3.121656 0.151589 0.032140 0.025873 0.002795 0.010250 0.070308 0.065669 0.016609 1.073790 0.040917 40.922701 15.426733 \
0.004063 0.079161 0.000000 0.112999 0.021444 0.371063 0.064924 2.075226 0.004177 0.037372 0.000155 0.004585 0.215788 0.007978 0.118229 0.016442 0.495176 10.291826 33.453780 15.127582 \
0.638696 0.001312 0.026551 0.040275 1.253945 0.002137 0.128111 0.073730 3.088481 0.340541 0.634065 0.001483 0.195073 0.664866 0.057328 0.438648 0.044742 0.775254 0.091276 0.286252 0.054021 \
0.000467 0.761771 0.000123 0.002163 0.000593 1.144692 0.014470 0.114551 0.265766 3.193996 0.000155 0.483076 0.369273 0.058614 0.617694 0.059927 0.330036 0.061583 0.730306 0.089835 0.364129 38.685701 \
0.126320 0.016628 0.576476 0.007508 0.508308 0.080383 2.066955 0.002179 0.486281 0.079236 0.163174 0.032232 0.055163 0.529045 0.071794 1.205738 0.033372 0.435109 0.074846 1.052040 0.063366 2.473439 0.751904 \
0.009760 0.107218 0.000000 0.250748 0.049246 0.423382 0.002122 1.519211 0.092070 0.332396 0.057910 0.105597 0.247490 0.079119 0.422671 0.105449 0.703795 0.107434 0.529594 0.184327 0.715716 1.106179 2.503268 17.923045 \
0.143832 0.000094 0.000741 0.003054 0.660622 0.001208 0.000579 0.001720 0.534375 0.001377 0.726908 0.077815 0.019696 0.663877 0.068758 0.134394 0.015019 0.500433 0.124232 0.063413 0.044676 2.460976 0.277265 1.164262 0.340811 \
0.000000 0.200806 0.000000 0.000064 0.000000 0.685812 0.000000 0.032106 0.000116 0.604541 0.012886 0.516927 0.176476 0.016022 0.544828 0.005436 0.563956 0.002398 0.563799 0.001702 0.798346 0.170088 2.478358 0.148940 2.029914 27.244097 \
0.030121 0.016020 0.136647 0.001527 0.006103 0.004089 0.557015 0.003211 0.043917 0.051686 0.232728 0.166150 0.146501 0.424607 0.112395 0.918198 0.041969 0.352807 0.154017 0.626603 0.091073 1.353860 0.526904 4.725840 0.617320 39.595443 12.677657 \
0.000934 0.027355 0.000000 0.127696 0.000085 0.004832 0.000000 1.903571 0.003713 0.081931 0.023909 0.183143 1.135910 0.039428 0.640495 0.040902 0.794366 0.009880 0.897101 0.010300 1.164525 0.316372 2.208430 0.299978 4.718199 12.868484 35.563093 30.574631 \
1.119411 0.059956 2.130663 1.292935 0.172403 0.000000 0.000386 0.000000 0.352731 0.000180 0.431456 0.000405 0.078312 3.330793 0.184010 1.328581 0.089308 0.292855 0.000096 0.002597 0.000246 0.193328 0.000000 0.078926 0.003859 0.076434 0.000000 0.000416 0.000000 \
0.056038 1.006045 0.042112 0.478019 0.000000 0.115975 0.000096 0.000344 0.000116 0.255975 0.000311 0.309643 0.136849 0.390190 3.765697 0.203017 2.469249 0.000096 0.270274 0.000448 0.021723 0.000469 0.127899 0.010543 0.105885 0.000118 0.238839 0.001248 0.003064 13.609310 \
1.075187 0.064968 5.159075 1.065537 0.000424 0.000000 0.403435 0.000000 0.573013 0.025454 0.069555 0.012138 0.170041 1.260239 0.136148 4.400610 0.048882 0.002014 0.000000 0.480521 0.000000 0.040109 0.000272 0.390087 0.000048 0.000000 0.000000 0.121855 0.000000 16.415611 5.784672 \
0.679370 0.800602 1.418466 3.062807 0.093491 0.042282 0.246094 0.527005 0.294368 0.300354 0.298091 0.324613 0.321642 1.220020 1.434579 1.635281 2.236557 0.081631 0.008455 0.042006 0.193459 0.323588 0.163406 0.443617 0.834976 0.028736 0.029786 0.015596 0.408680 1.155098 1.428293 2.230691 \
0.497293 0.000141 0.012473 0.015652 4.693944 0.487317 2.297807 0.199748 0.599932 0.000599 1.089585 0.001483 0.035939 0.831215 0.000060 0.004348 0.000070 1.050363 0.053805 0.345545 0.011476 0.898794 0.000272 0.374419 0.029088 0.344601 0.000000 0.001040 0.000000 1.266654 0.075878 0.351882 0.419831 \
0.000093 0.371541 0.000062 0.002990 0.196983 3.829580 0.150107 2.395833 0.000116 0.393545 0.000776 0.806071 0.037822 0.000396 0.897490 0.000679 0.073657 0.044125 0.870967 0.027138 0.527094 0.001125 0.969387 0.066519 0.617464 0.001060 1.481721 0.002079 0.017774 0.034573 1.066285 0.016701 0.433759 13.991583 \
0.079948 0.010539 0.871568 0.011134 2.018483 0.409721 5.709933 0.349846 0.208041 0.053363 0.415775 0.061227 0.060421 0.000857 0.000119 0.944560 0.000070 0.368538 0.026230 1.018005 0.015001 0.082373 0.021920 1.660885 0.001302 0.000589 0.000000 0.360575 0.000000 0.296014 0.048902 2.260424 0.853779 31.915858 8.373639 \
0.008032 0.036489 0.000062 0.586818 0.361164 1.562591 0.516594 4.919174 0.042119 0.177757 0.053874 0.173298 0.362681 0.000264 0.000595 0.000408 0.733202 0.079712 0.435243 0.083475 0.962295 0.076938 0.103080 0.002854 1.766961 0.004593 0.014243 0.002911 2.985421 0.090674 0.311759 0.154441 1.376727 12.116657 28.470047 19.459275 \
0.263567 0.000094 0.271628 0.077878 1.773102 0.000929 0.872084 0.040706 0.747870 0.042762 0.360038 0.000135 0.074859 0.259380 0.000000 0.019568 0.000140 0.787340 0.000192 0.096104 0.002705 2.691226 0.188587 1.759732 0.206851 0.682254 0.000068 0.009981 0.000981 0.239388 0.014351 0.256283 0.208924 2.057449 0.067554 1.243753 0.224397 \
0.000093 0.143614 0.002840 0.060699 0.003390 1.118208 0.187826 0.725836 0.046586 0.487814 0.000932 0.325422 0.053908 0.000330 0.187880 0.014540 0.034916 0.000959 0.532092 0.074161 0.095418 0.460970 2.203539 0.377447 0.985145 0.003180 0.863191 0.016636 0.065212 0.025405 0.175491 0.020837 0.219170 0.296066 1.346385 0.259909 0.822133 17.634677 \
0.148268 0.005996 0.612660 0.004963 0.340991 0.020909 1.496628 0.000459 0.467136 0.042942 0.099519 0.004316 0.051005 0.060922 0.002143 0.433620 0.000035 0.247195 0.012394 1.042457 0.000410 0.899262 0.143479 3.215817 0.285384 1.769879 0.296358 3.065510 0.011032 0.221536 0.023020 0.667397 0.275355 0.878163 0.089476 1.523251 0.199589 2.075154 0.413957 \
0.013122 0.043609 0.000062 0.333780 0.156468 0.251650 0.004438 0.768607 0.072867 0.140864 0.011489 0.034794 0.105226 0.039362 0.022384 0.000679 0.133032 0.220816 0.303613 0.012002 0.561523 0.324525 0.571469 0.461383 2.285052 0.560831 2.721043 0.034519 3.832200 0.041440 0.116405 0.056267 0.497593 0.291009 0.623366 0.256174 1.144639 0.524647 1.038682 12.931524 \
0.225554 0.000047 0.010929 0.009289 12.169045 0.083636 5.964323 0.681575 0.506470 0.000180 2.007768 0.181794 0.139046 0.322807 0.000060 0.009512 0.000105 1.035015 0.000288 0.021048 0.001066 1.293228 0.000453 1.177718 0.083261 0.772349 0.018146 0.530881 0.025006 0.359183 0.018727 0.269862 0.306375 4.943439 0.275865 2.397415 0.563566 4.971507 0.586685 1.293860 0.389004 \
0.000093 0.166706 0.000432 0.005790 0.094762 10.892976 1.049877 9.818281 0.000116 0.346890 0.248099 1.372357 0.167138 0.000330 0.300513 0.002718 0.017300 0.001247 1.337629 0.005195 0.104681 0.005060 1.282522 0.232701 1.418383 0.387941 1.320875 0.354545 1.360141 0.031140 0.238533 0.021539 0.304581 0.622868 4.699375 0.441084 2.871848 0.643789 4.127466 0.334224 0.928876 28.579806 \
0.140516 0.012600 0.423774 0.001718 0.047890 0.002044 1.094736 0.000115 0.424321 0.060909 0.144388 0.030883 0.145245 0.004747 0.000060 0.409432 0.000000 0.011511 0.000384 0.493508 0.000000 0.411115 0.025544 1.242140 0.000289 17.450524 1.113671 31.949764 2.418859 0.116039 0.012331 0.780008 0.305714 0.432918 0.021698 1.316696 0.087905 0.936840 0.273855 5.815294 1.197614 1.644621 0.403913 \
0.083310 0.056771 0.000247 0.506841 0.063994 0.028715 0.009261 0.514392 0.200499 0.269510 0.122186 0.070533 0.496706 0.009560 0.000952 0.000951 0.040320 0.060432 0.033244 0.009225 0.273876 0.140943 0.169022 0.003786 1.148387 6.798629 4.087042 15.287419 18.531553 0.093542 0.062369 0.317466 0.905810 0.309656 0.140701 0.511968 0.765495 0.454347 0.600415 1.868194 7.316623 1.477696 1.286990 43.916187 \
0.863970 0.065905 0.748196 0.529619 0.563995 0.000186 0.002219 0.000115 0.571505 0.000359 1.598824 0.004316 0.038763 1.897150 0.072866 0.555104 0.011580 0.708395 0.000192 0.009225 0.000246 0.338207 0.000091 0.150804 0.009600 0.336828 0.000000 0.003743 0.000000 6.546163 0.575921 2.577578 0.430124 2.898791 0.003020 0.056250 0.004868 0.318595 0.000295 0.201480 0.072138 0.501456 0.000219 0.032116 0.051709 \
0.026338 0.946136 0.005990 0.140867 0.000085 0.396897 0.000096 0.001261 0.000116 0.582801 0.001087 1.273908 0.092044 0.105361 2.720153 0.043756 1.085170 0.000192 0.760090 0.000537 0.031642 0.000375 0.491034 0.006757 0.069320 0.000589 0.648523 0.001871 0.005026 0.458622 7.487816 0.154207 0.350473 0.001027 2.862216 0.002515 0.010298 0.000000 0.215492 0.001930 0.056145 0.000097 0.381287 0.000205 0.002062 10.956917 \
0.565566 0.047403 2.299543 0.762425 0.001356 0.000279 0.813720 0.000115 0.236179 0.060430 0.754233 0.086986 0.058616 0.509595 0.031730 2.047159 0.020529 0.001151 0.000192 0.732111 0.000082 0.019586 0.002808 0.606945 0.000145 0.000353 0.000000 0.255147 0.000000 2.581654 0.313779 11.271062 0.569076 0.009008 0.001957 3.879355 0.002528 0.008844 0.002362 0.708204 0.000524 0.006305 0.000657 0.560642 0.002062 25.313949 5.637509 \
0.068927 0.371072 0.024699 1.108802 0.000254 0.001394 0.000772 0.451899 0.009630 0.116309 0.126844 0.530278 0.277072 0.091844 0.477439 0.173122 2.262490 0.000384 0.001537 0.000717 0.602428 0.032237 0.044112 0.000466 0.475496 0.001413 0.003287 0.000832 0.701399 0.953353 3.487342 0.911193 1.399673 0.003635 0.089643 0.011433 3.092500 0.000628 0.013582 0.000193 0.311885 0.001358 0.004160 0.000103 0.314379 8.832621 18.744445 13.945647 \
0.483563 0.000234 0.008892 0.035630 3.994417 0.771771 1.825878 0.250545 0.370309 0.000419 1.899865 0.011733 0.035860 0.879741 0.000417 0.004077 0.000772 1.272426 0.073021 0.346978 0.013526 0.420580 0.000091 0.272950 0.040907 0.324227 0.000000 0.001871 0.000000 0.536332 0.000253 0.001795 0.451134 2.359362 0.120121 0.324162 0.108501 0.643160 0.001329 0.216244 0.268662 2.323091 0.005328 0.013852 0.045374 2.600428 0.131929 0.662578 0.152215 \
0.000093 0.512203 0.000000 0.001654 0.111968 3.452570 0.084218 1.978448 0.000058 0.380788 0.000466 1.514097 0.044178 0.000132 1.191514 0.000136 0.103766 0.022062 1.394892 0.013077 0.971148 0.000562 0.696016 0.018814 0.634927 0.000236 1.610248 0.000416 0.020471 0.000040 0.536362 0.000000 0.169264 0.038559 2.159328 0.019665 0.498504 0.000045 0.565968 0.005500 0.373948 0.000485 2.214735 0.000000 0.001768 0.046583 2.620833 0.028569 0.682579 9.612709 \
0.109975 0.016535 1.041312 0.019406 1.931350 0.558500 4.380679 0.505677 0.176829 0.034737 0.806554 0.297371 0.031466 0.002374 0.000357 0.741135 0.000351 0.335924 0.069754 1.236457 0.087384 0.039265 0.004982 1.274118 0.002219 0.000589 0.000068 0.286547 0.000123 0.002262 0.000589 0.983926 0.517896 0.381796 0.077341 2.735831 0.318574 0.084620 0.041435 0.923644 0.004382 0.126382 0.063353 0.461729 0.004420 0.718719 0.092405 3.415722 0.415718 24.400553 6.746560 \
0.005884 0.074851 0.000000 0.220908 0.103323 1.262618 0.150589 4.658653 0.027035 0.106187 0.028567 0.586111 0.446015 0.000066 0.000893 0.000000 1.524024 0.014101 0.417565 0.017824 1.950083 0.080124 0.190037 0.001165 1.544626 0.001531 0.083744 0.000624 3.409178 0.000081 0.004629 0.000078 0.837302 0.023862 0.728891 0.049848 2.866325 0.003771 0.068501 0.000482 0.759132 0.006402 0.200205 0.000000 0.187832 0.054049 0.968351 0.081861 2.211488 5.140068 19.373137 11.561124 \
0.064397 0.000000 0.042112 0.038557 1.120532 0.003717 0.348448 0.117533 0.223763 0.015452 0.099985 0.000135 0.028249 0.129492 0.000000 0.012366 0.000491 0.661776 0.000769 0.147873 0.031560 0.746792 0.046739 0.706782 0.130873 0.162525 0.000000 0.007070 0.000368 0.066966 0.000042 0.001171 0.059065 0.928969 0.000559 0.092988 0.042595 3.529593 0.371685 0.604859 0.188097 1.702817 0.012481 0.030474 0.015763 0.153418 0.007112 0.078381 0.011491 0.396521 0.015140 0.189090 0.043198 \
0.000000 0.055366 0.000062 0.006808 0.000254 1.023142 0.007428 0.670108 0.010037 0.184704 0.000000 0.071612 0.066384 0.000066 0.135255 0.001359 0.015686 0.000096 0.976175 0.003672 0.644235 0.100928 0.975727 0.121389 0.928319 0.000236 0.915505 0.009981 0.150527 0.000000 0.032447 0.000000 0.011379 0.000158 1.013424 0.003354 0.095207 0.167041 2.729647 0.053168 0.426684 0.000388 2.005334 0.000718 0.008986 0.004101 0.119062 0.006776 0.041280 0.018617 0.802516 0.027912 0.702594 14.214694 \
0.084945 0.006464 0.287373 0.005472 0.330481 0.085680 1.265487 0.002179 0.257122 0.043721 0.028878 0.003641 0.009966 0.039560 0.002679 0.313495 0.000140 0.184749 0.105112 0.890822 0.005410 0.452442 0.106069 3.081614 0.536567 0.034978 0.025678 0.440217 0.000858 0.038612 0.009174 0.361403 0.033994 0.251423 0.109664 1.164866 0.003464 0.975582 0.193544 2.258321 0.308851 0.832592 0.308372 0.668173 0.004420 0.276499 0.042565 0.469281 0.055025 0.502355 0.140546 0.905488 0.227527 2.738552 0.892903 \
0.010974 0.034428 0.000000 0.159955 0.042380 0.283432 0.001061 1.029128 0.042815 0.136432 0.014439 0.013216 0.137634 0.004220 0.010061 0.000136 0.176300 0.034437 0.294294 0.001791 0.990330 0.159217 0.566034 0.343314 3.036767 0.007891 0.528692 0.001040 2.171984 0.003312 0.031984 0.000078 0.262465 0.033581 0.360196 0.000838 1.447392 0.149578 0.372719 0.159248 1.563846 0.129098 0.822643 0.000410 1.195790 0.049842 0.245019 0.053017 0.362328 0.106257 0.938586 0.157605 1.251589 1.091224 3.195698 12.984714 \
0.164659 0.000141 0.000741 0.003881 0.976185 0.001951 0.011673 0.007109 0.130940 0.000120 0.420899 0.045044 0.039313 0.169777 0.000060 0.000272 0.000175 0.418802 0.000288 0.002508 0.001312 0.388156 0.000091 0.042812 0.003377 0.241197 0.004656 0.042005 0.011768 0.069995 0.000000 0.000156 0.027479 0.380374 0.000112 0.000534 0.000374 1.322234 0.005905 0.048730 0.021649 2.382451 0.326035 0.037657 0.047437 0.164143 0.016776 0.072521 0.024883 1.572808 0.086923 0.585071 0.083552 0.629243 0.035120 0.089148 0.030223 \
0.000000 0.172889 0.000000 0.000191 0.000085 0.880032 0.000289 0.356038 0.000058 0.127388 0.007608 0.309374 0.105305 0.000000 0.240505 0.000000 0.047268 0.000096 0.636916 0.000090 0.395771 0.000843 0.566759 0.016193 0.336277 0.021435 0.676049 0.008942 0.703728 0.000283 0.055425 0.000000 0.018603 0.000000 0.518903 0.000000 0.006459 0.001122 1.110726 0.002863 0.176224 0.054025 2.392606 0.000821 0.012227 0.002050 0.201477 0.001557 0.051048 0.022214 1.797671 0.027973 1.398079 0.037461 1.228004 0.030585 0.239725 13.935950";
string model_ECMunrest = model_ECMunrest1 + 
"0.113991 0.018315 0.201112 0.001082 0.012121 0.001951 1.720919 0.001720 0.082323 0.029826 0.197641 0.061497 0.073682 0.000330 0.000060 0.165784 0.000070 0.003549 0.000384 0.556204 0.000164 0.097554 0.004982 0.551493 0.000289 0.015310 0.000753 0.247245 0.010419 0.000283 0.000084 0.194319 0.037724 0.002449 0.000112 0.466770 0.000187 0.909861 0.280400 0.713961 0.001760 1.179053 0.298738 0.938439 0.165587 0.080337 0.009773 0.324696 0.016839 0.658541 0.036022 1.693998 0.046588 0.375097 0.067431 0.639125 0.053748 21.171295 5.214689 \
0.018773 0.032039 0.000000 0.175861 0.002797 0.002974 0.003376 2.163175 0.007948 0.014314 0.105884 0.183952 0.381671 0.000066 0.000119 0.000000 0.185038 0.001918 0.001441 0.001254 0.703092 0.084060 0.053714 0.003029 0.634203 0.043222 0.097165 0.143481 0.590833 0.000081 0.000295 0.000078 0.410199 0.000553 0.000447 0.000610 0.716441 0.194964 0.293884 0.001158 0.744000 0.684968 1.149846 0.069567 1.558784 0.032177 0.064227 0.074536 0.276276 0.238907 0.496552 0.672077 1.526141 0.235747 0.403521 0.136937 0.968146 13.981617 18.675227 25.640860 \
\
0.021414 0.021349 0.016195 0.015717 0.011798 0.010761 0.010366 0.008721 0.017237 0.016697 0.006441 0.007415 0.012744 0.015167 0.016798 0.007359 0.028497 0.010425 0.010408 0.011165 0.012199 0.010671 0.011040 0.017168 0.020730 0.008491 0.014604 0.004809 0.008158 0.024759 0.023762 0.012814 0.021180 0.012656 0.017882 0.013120 0.010682 0.022276 0.020321 0.031090 0.026699 0.010310 0.013701 0.009746 0.006788 0.019020 0.018419 0.010921 0.022626 0.018907 0.026817 0.016516 0.018288 0.028590 0.025285 0.034527 0.030606 0.016883 0.023659 0.016386 0.010223 \
\
TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
GGG";

/* empirical codon model of Schneider, Cannarozzi, Gonnet 2005
*/

string model_ECM_Schneider05 = 
  "15594\
    787    609\
    717    391  22864\
    476      0     89    133\
      0    444     39     51  15656\
     77     11   2230    662  15154  11009\
     99      0      0   2501  21330  19335  36566\
   2355    360    127     41    400     18     28      0\
    300   1911     76     49      0    260     23     70  25720\
    378     30     81     96   1004      0    238     90   1025      0\
     60    275     46     53     37    766     83    157      0    623  29318\
    102     74     24    357     16     23     83    277    110     80    186    164\
   1304      0   4065   4641    278      0      0      0    138     18    147     21      9\
      0   1293   3595   2691      0    214      0      0      1    128      9    152      9  15530\
    112    106  26307   6745      0      0    154    340     68      0     59      0     36  10903   8429\
    186    136   4180  13314     20     34      0      0     39     14     17      9     51   7205   7138  16270\
     97      0      7     72   2744      0      0      0     43      3    101     12      0    937     34      0      0\
      0     82     18      3     65   2535      0      0      0     47      0     49      0      0    706      0      0  16327\
      0     12     99     87      0      0   2829    131      7      0     11      0     14      0      0   1474    137  15604  10490\
     26      3     11    109      0      0    244   5851      4     22     44     31     68      0      0      0    935  15588  21512  27394\
    172     22     80     27    315     22     76      0   2184      0    475      6     63    676     55      0     43    756     38     56    173\
     36    166      5     20     43    211     53    112    137   1624     21    276     20      0    390     96     40      0    617     97    246  30148\
     35     18     93     63    143     90    403      0    122     56     70     60     38    168     38   1166     51    126     99   1405    296   1995   1176\
     11     13     50     95     55     80    109    288     42     58     43     35     97     77     92     46    290     88    161    235   1234   1086   1057  15317\
     70      7     49     98    326     11      1      0    362      0   1302      0     14    685      0      0      0    712      0      0      0   4390      0    646    232\
      0     56     40      0      0    177      0    105      4    139      0    908     40      1    486      0     34      0    404      8    121      0   2689    217    355  42123\
      8      0      0     25      0      0    163     53     68      0     51     53    130      0     45    966      0      0     23    393      0    133    261   3956      0  40836  19396\
     10      8      2     45      0     50     31    385      0     17     50    109    576     21      0      7    340      0      0     57   1276    191    317      0   2408  25507  19595  39593\
    437     29    201    399     72     30      0     18     64      1     90      1     10   1848    131    185    283     68     10     24      0     84      4      0     20     27      0     64      0\
     23    356    111     27     15     35      6     15      7     42     10     62      7    133   1745      0    351      0     55      0     29      0     36     24      7     17     20      0     15  16920\
    157     33   3243    833      0     54    209     88     29      1     14     16     14    134    176   3786    421     49     13     97     79     51      4     93     20     43      3      5      0  15113  13126\
    126     89    536   1607     75     39    153    229     32     21     29     22     45    346    369    726   1104     34     27     55     92     61     20     96    116     29     23      0     50    781    495   2443\
     82     23     67     79   2679    170    489    362     61     17    194     41      0    233      0      0     25    990     52      0     15    212      6    166     30    318     42      0      0   1552      0      0    262\
      7     74      3     39    279   2315    113      0     21     33     10    153      9      0    186      0     44     76    817      0      0     15    171     37     61     14    118     71     33      0   1028      0    176  16316\
     23     11    295     92    567     81   3129    578     12      0     75     22     19      4      0    288     46     64      0   1015     88     93     27    341     96     39      0      0     26      0      0   2434    714  16829  10306\
     51      5     49    306      0    682    232   6944     52      9    134      0     44     55     28      0    181      0    170      0   1953      0     35      0    307      0    116    322    172    166    134      0   1616  20814  19942  34110\
     41      0     48     36    289     32    139     78    316      0    191      0     11    110      0     22      0    115     34     55     46   1623     55    311    168    429      0     45     27    202      0     50     68   1425     16    237    182\
      8     34      4     22     87    277     79     40      0    196      4    116      7      0     64     19      7     14     98     17     47      9   1196    148    142      0    264      0     58      0    145     21     42    122    984    165     67  17900\
      4      6     59     28     54     24    133     53     16      8     20     15      9     30     17    100      0     26     18     91     39    169     94   1364    137    334    337    908    410     16      0    292     80    114     54    650      0    580    305\
      6      0      4     39     26     24     52     70     11     11     13     16     24      7      3      0     36     14     14     28    105     84    110     33    727    353    569    299   1249     20     21     22    189     27     85     94    795    298    269  12057\
     49      0      0     47    938     51    528    366    157      0    981    105     27     71      0     26      4    178      0     76     37    605      0    198     84   1145      0      0     44    345      0     91     78   2941      0    402    387   4830      0    245    146\
     10     34     33      0    312    983    214    697     10     79    169    792     19     23     68      0     21     15    213     25     88     41    454    133    123      0    740     72     22      0    228     27     67    421   2248    343    618    178   3430    158    131  21641\
     10      2     87     13     54     26     79     12      0     20     52      7     23     32      0      9     13      5     17      0     45    118     55    623     32   6203   2803  22415   6405      0      0    415     93     11     31    743      0    178    142   4144    390    638    339\
     12      0     51     36     33     12      4     72     31      0     36     20    264      0     20      0     22     17      0     39      5     62     95     60    271   2885   5913   6064  18118     29     15    106    258     40     59     49    927    145    139    572   3644    544    444  25118\
    446      0    175    236    218     61     96      0     80      0    108     60     20   1661     28      0     11    258      0      0     65     68      1     82     28     96     15      0      3   6623     53    389    259   1227     72     65    141    100     20     37      3    233     69     31      0\
      6    460     73     35     15    199     49      0      7     69     28    105      4     18   1603    120     42      0    199     36     69     15     75     29     25     18     52      0     36     79   6084    158    190    148   1080      1      0     21     86     13     11     16    170      5     38  17102\
     83     13   3202    251     70      2    396      0     18      2     52     24      6    131      5   2164      7      0     63    394      0     23     36    276     31      0      0     86      0    125      0  11512    626      0     20   1694     36     55     26    153     16     31     87    249      0  15228  10722\
     60     42      0   1017     46     19     85    271      9     17     54     16     27      0     24     16    827     27     21     58    154     23     18     10     62     12     20      0     58    743    928   1567   1219    155     66    206    995     12     12     10     54     28     10      0    106  10663   9288  21679\
     89     10     43     68   2454    180    301      0     54     11    203     32      0    183      8     14     31   1111     63     28      0    111      5     60     52    121      0      0     62    233     32     42    108   3490      0      0      0    252     38     43     31    628    101     47      0   2567      0     14     37\
      1     78     42     24    250   2338    201    324     11     28     22    197     15     30    160      0     31     79    965      7    145      0    119     22     88     29    133     26      0      0    217      0    115      0   2888      0     36     25    224     13     45      0    608     24     51      0   2241      0     80  12825\
     39      0    293     82    303    127   3185    250     19      1     87     27     17     29      0    300     46     91     19   1186    101     70     20    307     61     38      0     72     23     73      0    484    282      0      0   3718      0    149     28    216     48    268    106    108     88      0      0   3240    591  14816   8169\
     19     59      0    480    450    354      0   9766      0     68     96    134     32     16     40      0    261      0    206      0   4882    109     45    182    246      0     94      0    379      0     77    315    392      0    193      0   7718     81    161     62    171    178    258     49    171      0      0      0   2849  13866  14245  22755\
     14      0      0     16    138      0     52     20    122      0     41      0      3     21      0      5      0     44      7     28     37    378      0    132     54     66      8      0     13     36      0     19      8    149      2     83      0   1636      4     56     44    488     30     27     19    187      0     24      5    384      0     59     20\
      0     13     16      0      1    115     35     76      0     76      9     24      1     10     20      0      3      8     72      6     40     22    311     76     72      1     57      0     11      0     20      0      9     54    124      0     77     33   1249     58     35     38    395     22     11      0    141     19      0      0    322     39    182  13482\
      8      2     32     46     42     14    118     53     16      9     13      0      5     13      6     52      0     29     17     71     35    123     28   1246    112     69      0     24      0      0      5     95     22     38     36    196      0    187     74    626      0    114     82    213      0     62     24    367     16     64     41    636      0   1758   1114\
      3      4     20     18     36     42     43     84     16      6      3     11     13      7     12      0     27     19     40     51     81     69     71     66    739      0     44      0    104     21      1     17     63     52     35     76    199    106     80      0    429    101     86      0    162      4     26     38    161     70     91    118    916   1216   1466  11181\
     32      0     44      0    293     18     60     85     62      0    515      0      0     50      0      9      9    161     16     29      0    266      0      0     43    688      0    117      0     46      0      0     23    349     91     46      0    829     98     44     32   3432      6      0      0    557      0      0     21   1264      0    171      0   1058      0     28     30\
      9     23      0     15     52    224     28    240      0     38     38    372      7     14     38      0      6     14    114     31    129     12    165     15     40      1    452      0     78     10     42     18     16     50    250      6     62     91    589     13     28      0   2627     25      0     15    383      0      0     72   1050      0    357      0    795      0     50  17792\
     15      1     23      0     80     59    198      0      9      0     34      0      4     13     12     70      0     16     11     78     26     26     16    225      8      0      0    176      0     16      0     85     31      3     12    214     51    121     31    137     11    151    134   1071      0      0      0    629      0    133     44   1071      0     50     86    798      0  13656   9267\
      2     11      9    122     10     49    162    532     14      3     42     39    214      0      1      0     39      9     55     36    251     29     31     53    148      0      0      0    515      0     22     14    116     77      0     84    405     67    171     13    144    398    328      0   1567      0      0     86    380     35    151    360   2466     88     32      0    749  11092  11188  21307\
 0.019065 0.019572 0.008521 0.014125 0.017115 0.015966 0.013340 0.004325 0.013196 0.015988 0.010982 0.011900 0.012462 0.015017 0.017340 0.008031 0.037328 0.017599 0.014989 0.017755 0.005879 0.011592 0.014372 0.013669 0.033679 0.005419 0.008712 0.006185 0.008488 0.017469 0.019956 0.009364 0.021871 0.014391 0.015964 0.016870 0.005890 0.018236 0.020614 0.028227 0.031883 0.013622 0.019058 0.013502 0.011845 0.013592 0.013760 0.008372 0.026506 0.019952 0.021260 0.017896 0.005964 0.025167 0.024586 0.031093 0.038868 0.011549 0.017608 0.018437 0.013269\
TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
GGG";


ModelCodon::ModelCodon(const char *model_name, string model_params, StateFreqType freq, string freq_params,
		PhyloTree *tree) : ModelMarkov(tree)
{
    half_matrix = false;
    omega = kappa = kappa2 = 1.0;
    fix_omega = fix_kappa = false;
    fix_kappa2 = true;
    codon_freq_style = CF_TARGET_CODON;
    codon_kappa_style = CK_ONE_KAPPA;
	ntfreq = new double[12];
	empirical_rates = NULL;
	int nrates = getNumRateEntries();
    delete [] rates;
    rates = new double[nrates];
    empirical_rates = new double [nrates];

    rate_attr = NULL;
    computeRateAttributes();

   	init(model_name, model_params, freq, freq_params);
}

ModelCodon::~ModelCodon() {
	if (rate_attr) {
		delete [] rate_attr;
		rate_attr = NULL;
	}
	if (empirical_rates) {
		delete [] empirical_rates;
		empirical_rates = NULL;
	}
	if (ntfreq) {
		delete [] ntfreq;
		ntfreq = NULL;
	}
}

void ModelCodon::startCheckpoint() {
    checkpoint->startStruct("ModelCodon");
}

void ModelCodon::saveCheckpoint() {
    startCheckpoint();
//    CKP_ARRAY_SAVE(12, ntfreq);
    CKP_SAVE(omega);
//    CKP_SAVE(fix_omega);
//    int codon_kappa_style = this->codon_kappa_style;
//    CKP_SAVE(codon_kappa_style);
    CKP_SAVE(kappa);
//    CKP_SAVE(fix_kappa);
    CKP_SAVE(kappa2);
//    CKP_SAVE(fix_kappa2);
//    int codon_freq_style = this->codon_freq_style;
//    CKP_SAVE(codon_freq_style);
    endCheckpoint();

    ModelMarkov::saveCheckpoint();
}

void ModelCodon::restoreCheckpoint() {
    ModelMarkov::restoreCheckpoint();
    startCheckpoint();
//    CKP_ARRAY_RESTORE(12, ntfreq);
    CKP_RESTORE(omega);
//    CKP_RESTORE(fix_omega);
//    int codon_kappa_style;
//    CKP_RESTORE(codon_kappa_style);
//    this->codon_kappa_style = (CodonKappaStyle)codon_kappa_style;
    CKP_RESTORE(kappa);
//    CKP_RESTORE(fix_kappa);
    CKP_RESTORE(kappa2);
//    CKP_RESTORE(fix_kappa2);
//    int codon_freq_style;
//    CKP_RESTORE(codon_freq_style);
//    this->codon_freq_style = (CodonFreqStyle)codon_freq_style;
    endCheckpoint();

    decomposeRateMatrix();
    if (phylo_tree)
        phylo_tree->clearAllPartialLH();

}

StateFreqType ModelCodon::initCodon(const char *model_name, StateFreqType freq, bool reset_params) {
	string name_upper = model_name;
	for (string::iterator it = name_upper.begin(); it != name_upper.end(); it++)
		(*it) = toupper(*it);
    
	if (name_upper == "MG") {
		return initMG94(true, freq, CK_ONE_KAPPA);
	} else if (name_upper == "MGK") {
		return initMG94(false, freq, CK_ONE_KAPPA);
	} else if (name_upper == "MG1KTS" || name_upper == "MGKAP2") {
        return initMG94(false, freq, CK_ONE_KAPPA_TS);
	} else if (name_upper == "MG1KTV" || name_upper == "MGKAP3") {
        return initMG94(false, freq, CK_ONE_KAPPA_TV);
	} else if (name_upper == "MG2K" || name_upper == "MGKAP4") {
        return initMG94(false, freq, CK_TWO_KAPPA);
	} else if (name_upper == "GY") {
        return initGY94(false, CK_ONE_KAPPA);
	} else if (name_upper == "GY0K" || name_upper == "GYKAP1") {
        return initGY94(true, CK_ONE_KAPPA);
	} else if (name_upper == "GY1KTS" || name_upper == "GYKAP2") {
        return initGY94(false, CK_ONE_KAPPA_TS);
	} else if (name_upper == "GY1KTV" || name_upper == "GYKAP3") {
        return initGY94(false, CK_ONE_KAPPA_TV);
	} else if (name_upper == "GY2K" || name_upper == "GYKAP4") {
        return initGY94(false, CK_TWO_KAPPA);
	} else if (name_upper == "ECM" || name_upper == "KOSI07" || name_upper == "ECMK07") {
		if (!phylo_tree->aln->isStandardGeneticCode())
			outError("For ECMK07 a standard genetic code must be used");
		readCodonModel(model_ECMunrest, reset_params);
		return FREQ_USER_DEFINED;
	} else if (name_upper == "ECMREST") {
		if (!phylo_tree->aln->isStandardGeneticCode())
			outError("For ECMREST a standard genetic code must be used");
		readCodonModel(model_ECMrest, reset_params);
		return FREQ_USER_DEFINED;
	} else if (name_upper == "SCHN05" || name_upper == "ECMS05") {
		if (!phylo_tree->aln->isStandardGeneticCode())
			outError("For ECMS05 a standard genetic code must be used");
		readCodonModel(model_ECM_Schneider05, reset_params);
		return FREQ_USER_DEFINED;
	} else {
		//cout << "User-specified model "<< model_name << endl;
		//readParameters(model_name);
			//name += " (user-defined)";
        readCodonModelFile(model_name, reset_params);
		return FREQ_USER_DEFINED;
	}

	return FREQ_UNKNOWN;
}

void ModelCodon::init(const char *model_name, string model_params, StateFreqType freq, string freq_params)
{
    int i, j;
	for (i = 0; i < 12; i++)
		ntfreq[i] = 0.25;
    // initialize empirical_rates
    for (i = 0; i < num_states; i++) {
        double *this_emp_rate = &empirical_rates[i*num_states];
        int *this_rate_attr = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            memset(this_emp_rate, 0, num_states*sizeof(double));
            continue;
        }
        for (j = 0; j < num_states; j++) {
            int attr = this_rate_attr[j];
            if (attr & (CA_STOP_CODON+CA_MULTI_NT)) { // stop codon or multiple nt substitutions
                this_emp_rate[j] = 0.0;
            } else {
                this_emp_rate[j] = 1.0;
            }
        }
    }    

    ignore_state_freq = false;

	StateFreqType def_freq = FREQ_UNKNOWN;
	name = full_name = model_name;
    size_t pos;
	if ((pos=name.find('_')) == string::npos) {
		def_freq = initCodon(model_name, freq, true);
	} else {
		def_freq = initCodon(name.substr(0, pos).c_str(), freq, false);
		if (def_freq != FREQ_USER_DEFINED)
			outError("Invalid model " + name + ": first component must be an empirical model"); // first model must be empirical
		def_freq = initCodon(name.substr(pos+1).c_str(), freq, false);
		if (def_freq == FREQ_USER_DEFINED) // second model must be parametric
			outError("Invalid model " + name + ": second component must be a mechanistic model");
		// adjust the constraint
        if (codon_freq_style==CF_TARGET_CODON) 
            def_freq = FREQ_USER_DEFINED;
	}

    num_params = (!fix_omega) + (!fix_kappa) + (!fix_kappa2);

	if (freq_params != "") {
		readStateFreq(freq_params);
	}
	if (model_params != "") {
		readRates(model_params);
	}

//	if (freq == FREQ_UNKNOWN ||  def_freq == FREQ_EQUAL) freq = def_freq;
    if (freq == FREQ_UNKNOWN) freq = def_freq;
	if (freq == FREQ_CODON_1x4 || freq == FREQ_CODON_3x4 || freq == FREQ_CODON_3x4C) {
		//ntfreq = new double[12];
		phylo_tree->aln->computeCodonFreq(freq, state_freq, ntfreq);
	}
	ModelMarkov::init(freq);
}

StateFreqType ModelCodon::initMG94(bool fix_kappa, StateFreqType freq, CodonKappaStyle kappa_style) {
	/* Muse-Gaut 1994 model with 1 parameters: omega */

    fix_omega = false;
    this->fix_kappa = fix_kappa;
    if (fix_kappa)
        kappa = 1.0;
    fix_kappa2 = true;
    codon_freq_style = CF_TARGET_NT;
    this->codon_kappa_style = kappa_style;
    if (kappa_style == CK_TWO_KAPPA)
        fix_kappa2 = false;
    
    if (freq == FREQ_UNKNOWN || freq == FREQ_USER_DEFINED)
        freq = FREQ_CODON_3x4;
        
    switch (freq) {
      case FREQ_CODON_1x4:
      case FREQ_CODON_3x4:
      case FREQ_CODON_3x4C:
		phylo_tree->aln->computeCodonFreq(freq, state_freq, ntfreq);
        break;
      case FREQ_EMPIRICAL:
      case FREQ_ESTIMATE:
      case FREQ_USER_DEFINED:
        outError("Invalid state frequency type for MG model, please use +F1X4 or +F3X4 or +F3X4C");
        break;
      default:
        break;
    }
    
    // ignote state_freq because ntfreq is already used
    ignore_state_freq = true;
    combineRateNTFreq();
    
    return FREQ_CODON_3x4;
}

StateFreqType ModelCodon::initGY94(bool fix_kappa, CodonKappaStyle kappa_style) {
    fix_omega = false;
    this->fix_kappa = fix_kappa;
    if (fix_kappa)
        kappa = 1.0;
    fix_kappa2 = true;
    codon_freq_style = CF_TARGET_CODON;
    this->codon_kappa_style = kappa_style;
    if (kappa_style == CK_TWO_KAPPA)
        fix_kappa2 = false;
            
    return FREQ_EMPIRICAL;
}


void ModelCodon::computeRateAttributes() {
    char symbols_protein[] = "ARNDCQEGHILKMFPSTWYVX"; // X for unknown AA
    int i, j, ts, tv;
    int nrates = getNumRateEntries();
    if (!rate_attr) {
        rate_attr = new int[nrates];
        memset(rate_attr, 0, sizeof(int)*nrates);
    }
    char aa_cost_change[400];
    memset(aa_cost_change, 10, sizeof(char)*400);
    for (i = 0; i < num_states; i++) {
        int *rate_attr_row = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            for (j = 0; j < num_states; j++)
                rate_attr_row[j] = CA_STOP_CODON;
            continue;
        }
        for (j = 0; j < num_states; j++)  {
            if (j == i || phylo_tree->aln->isStopCodon(j)) {
                rate_attr_row[j] = CA_STOP_CODON;
                continue;
            }
            int nuc1, nuc2;
            int attr = 0;
            ts = tv = 0;
            int codoni = phylo_tree->aln->codon_table[i];
            int codonj = phylo_tree->aln->codon_table[j];
            if (phylo_tree->aln->genetic_code[codoni] == phylo_tree->aln->genetic_code[codonj])
                attr |= CA_SYNONYMOUS;
            else
                attr |= CA_NONSYNONYMOUS;
                
                
            int nt_changes = ((codoni/16) != (codonj/16)) + (((codoni%16)/4) != ((codonj%16)/4)) + ((codoni%4) != (codonj%4));
            int aa1 = strchr(symbols_protein, phylo_tree->aln->genetic_code[codoni]) - symbols_protein;
            int aa2 = strchr(symbols_protein, phylo_tree->aln->genetic_code[codonj]) - symbols_protein;
            ASSERT(aa1 >= 0 && aa1 < 20 && aa2 >= 0 && aa2 < 20);
            if (nt_changes < aa_cost_change[aa1*20+aa2]) {
                aa_cost_change[aa1*20+aa2] = aa_cost_change[aa2*20+aa1] = nt_changes;
            }
            
            if ((nuc1=codoni/16) != (nuc2=codonj/16)) {
                if (abs(nuc1-nuc2)==2) { // transition 
                    attr |= CA_TRANSITION_1NT;
                    ts++;
                } else { // transversion
                    attr |= CA_TRANSVERSION_1NT;
                    tv++;
                }
            }
            if ((nuc1=(codoni%16)/4) != (nuc2=(codonj%16)/4)) {
                if (abs(nuc1-nuc2)==2) { // transition
                    attr |= CA_TRANSITION_2NT;
                    ts++;
                } else { // transversion
                    attr |= CA_TRANSVERSION_2NT;
                    tv++;
                }
            }
            if ((nuc1=codoni%4) != (nuc2=codonj%4)) {
                if (abs(nuc1-nuc2)==2) { // transition
                    attr |= CA_TRANSITION_3NT;
                    ts++;
                } else { // transversion
                    attr |= CA_TRANSVERSION_3NT;
                    tv++;
                }
            }
            if (ts+tv>1) 
                attr |= CA_MULTI_NT;
            else if (ts==1) 
                attr |= CA_TRANSITION;
            else if (tv==1)
                attr |= CA_TRANSVERSION;
                    
            rate_attr_row[j] = attr;
        }
    }
    
    if (verbose_mode >= VB_MAX) {

        // make cost matrix fulfill triangular inequality
        for (int k = 0; k < 20; k++)
            for (i = 0; i < 20; i++)
                for (j = 0; j < 20; j++)
                    if (aa_cost_change[i*20+j] > aa_cost_change[i*20+k] + aa_cost_change[k*20+j])
                        aa_cost_change[i*20+j] = aa_cost_change[i*20+k] + aa_cost_change[k*20+j];

        cout << "cost matrix by number of nt changes for TNT use" << endl;
        cout << "smatrix =1 (aa_nt_changes)";
        for (i = 0; i < 19; i++)
            for (j = i+1; j < 20; j++)
                cout << " " << symbols_protein[i] << "/" << symbols_protein[j] << " " << (int)aa_cost_change[i*20+j];
        cout << ";" << endl;
        cout << 20 << endl;
        for (i = 0; i < 20; i++) {
            aa_cost_change[i*20+i] = 0;
            for (j = 0; j < 20; j++)
                cout << (int)aa_cost_change[i*20+j] << " ";
            cout << endl;
        }
        
    }
    
}

void ModelCodon::combineRateNTFreq() {
    int i, j;
    for (i = 0; i < num_states; i++) {
        if (phylo_tree->aln->isStopCodon(i))
            continue;
        double *this_rate = &empirical_rates[i*num_states];
        for (j = 0; j < num_states; j++)  {
            if (this_rate[j] == 0.0)
                continue;

            int codoni = phylo_tree->aln->codon_table[i];
            int codonj = phylo_tree->aln->codon_table[j];
            int nuc1, nuc2;
                
            if ((nuc1=codoni/16) != (nuc2=codonj/16)) {
                this_rate[j] *= ntfreq[nuc2];
            }
            if ((nuc1=(codoni%16)/4) != (nuc2=(codonj%16)/4)) {
                this_rate[j] *= ntfreq[nuc2+4];
            }
            if ((nuc1=codoni%4) != (nuc2=codonj%4)) {
                this_rate[j] *= ntfreq[nuc2+8];
            }
        }
    }
    
}


void ModelCodon::readCodonModel(istream &in, bool reset_params) {
	int nrates = getNumRateEntries();

	int i, j;
	int nscodons = phylo_tree->aln->getNumNonstopCodons();

	double * q = new double[nscodons*nscodons];
	double *f = new double[nscodons];
	for (i = 1; i < nscodons; i++) {
		for (j = 0; j < i; j++) {
			in >> q[i*nscodons+j];
			//q[j*num_states+i] = q[i*num_states+j];
			if (verbose_mode >= VB_MAX) cout << " " << q[i*nscodons+j];
		}
		if (verbose_mode >= VB_MAX) cout << endl;
	}
	for (i = 0; i < nscodons; i++)
		in >> f[i];
	StrVector codons;
	codons.resize(nscodons);
	IntVector state_map;
	state_map.resize(nscodons);
	for (i = 0; i < nscodons; i++) {
		in >> codons[i];
		if (codons[i].length() != 3)
			outError("Input model has wrong codon format ", codons[i]);
		int nt1 = phylo_tree->aln->convertState(codons[i][0], SEQ_DNA);
		int nt2 = phylo_tree->aln->convertState(codons[i][1], SEQ_DNA);
		int nt3 = phylo_tree->aln->convertState(codons[i][2], SEQ_DNA);
		if (nt1 > 3 || nt2 > 3 || nt3 > 3)
			outError("Wrong codon triplet ", codons[i]);
		state_map[i] = phylo_tree->aln->non_stop_codon[nt1*16+nt2*4+nt3];
		if (phylo_tree->aln->isStopCodon(state_map[i]) || state_map[i] == STATE_INVALID)
			outError("Stop codon encountered");
		if (verbose_mode >= VB_MAX)
			cout << " " << codons[i] << " " << state_map[i];
	}
	if (verbose_mode >= VB_MAX) cout << endl;

	//int row = 0, col = 1;
	// since rates for codons is stored in lower-triangle, special treatment is needed
    memset(empirical_rates, 0, nrates*sizeof(double));
    memset(rates, 0, nrates*sizeof(double));
	for (i = 1; i < nscodons; i++) {
		for (j = 0; j < i; j++) {
			int row = state_map[i], col = state_map[j];
			if (row < col) {
				int tmp = row;
				row = col;
				col = tmp;
			}
//			int id = col*(2*num_states-col-1)/2 + (row-col-1);
            double qentry = q[i*nscodons+j];
            int id = row*num_states+col;
			ASSERT(id < nrates && id >= 0);
			empirical_rates[id] = rates[id] = qentry;
            id = col*num_states+row;
			ASSERT(id < nrates && id >= 0);
			empirical_rates[id] = rates[id] = qentry;
		}
	}
	memset(state_freq, 0, num_states*sizeof(double));
	for (i = 0; i < num_states; i++)
		state_freq[i] = MIN_FREQUENCY;
	for (i = 0; i < nscodons; i++)
		state_freq[state_map[i]] = f[i]-(num_states-nscodons)*MIN_FREQUENCY/nscodons;

    if (reset_params) {
        fix_omega = fix_kappa = fix_kappa2 = true;
        omega = kappa = kappa2 = 1.0;
    }
	delete [] f;
	delete [] q;
}

void ModelCodon::readCodonModel(string &str, bool reset_params) {
	try {
		istringstream in(str);
		readCodonModel(in, reset_params);
	}
	catch (const char *str) {
		outError(str);
	}
}

void ModelCodon::readCodonModelFile(const char *filename, bool reset_params) {
	try {
		ifstream in;
        // set the failbit and badbit
        in.exceptions(ios::failbit | ios::badbit);
        in.open(filename);
        // remove the failbit
        in.exceptions(ios::badbit);

		readCodonModel(in, reset_params);

        in.clear();
        // set the failbit again
        in.exceptions(ios::failbit | ios::badbit);
        in.close();
	}
	catch (const char *str) {
		outError(str);
	}
	catch (...) {
        outError(ERR_READ_INPUT, filename);
    }
}

void ModelCodon::decomposeRateMatrix() {
    computeCodonRateMatrix();
    ModelMarkov::decomposeRateMatrix();
}

void ModelCodon::computeCodonRateMatrix() {
//    if (num_params == 0) 
//        return; // do nothing for empirical codon model
        
    switch (codon_kappa_style) {
    case CK_ONE_KAPPA:
        computeCodonRateMatrix_1KAPPA();
        break;
    case CK_ONE_KAPPA_TS:
        computeCodonRateMatrix_1KAPPATS();
        break;
    case CK_ONE_KAPPA_TV:
        computeCodonRateMatrix_1KAPPATV();
        break;
    case CK_TWO_KAPPA:
        computeCodonRateMatrix_2KAPPA();
        break;
    }
}

void ModelCodon::computeCodonRateMatrix_1KAPPA() {
    int nrates = getNumRateEntries();
    memcpy(rates, empirical_rates, nrates*sizeof(double));
    if (omega == 1.0 && kappa == 1.0)
        return; // do nothing

    int i, j;
    double omega_kappa = omega*kappa;
    
    for (i = 0; i < num_states; i++) {
        double *this_rate = &rates[i*num_states];
        int *this_rate_attr = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            continue;
        }
        for (j = 0; j < num_states; j++) {
            if (this_rate[j] == 0.0) continue;
            int attr = this_rate_attr[j];
            if (attr & CA_SYNONYMOUS) { // synonymous
                if (attr & CA_TRANSITION) // transition
                    this_rate[j] *= kappa;
            } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
                if (attr & CA_TRANSITION) // transition
                    this_rate[j] *= omega_kappa;                
                else // transversion
                    this_rate[j] *= omega;
            }
        }
    }
}

void ModelCodon::computeCodonRateMatrix_1KAPPATS() {
    int nrates = getNumRateEntries();
    memcpy(rates, empirical_rates, nrates*sizeof(double));

    int i, j;
    double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
    double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};

    for (i = 0; i < num_states; i++) {
        double *this_rate = &rates[i*num_states];
        int *this_rate_attr = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            continue;
        }
        for (j = 0; j < num_states; j++) {
            int attr = this_rate_attr[j];
            if (this_rate[j] == 0.0) continue;
            if (attr & CA_SYNONYMOUS) { // synonymous
                int num = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
                this_rate[j] *= kappa_pow[num];
            } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
                int num = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
                this_rate[j] *= omega_kappa_pow[num];
            }
        }
    }
}

void ModelCodon::computeCodonRateMatrix_1KAPPATV() {
    int nrates = getNumRateEntries();
    memcpy(rates, empirical_rates, nrates*sizeof(double));

    int i, j;
    double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
    double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};

    for (i = 0; i < num_states; i++) {
        double *this_rate = &rates[i*num_states];
        int *this_rate_attr = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            continue;
        }
        for (j = 0; j < num_states; j++) {
            int attr = this_rate_attr[j];
            if (this_rate[j] == 0.0) continue;
            if (attr & CA_SYNONYMOUS) { // synonymous
                int num = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
                this_rate[j] *= kappa_pow[num];
            } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
                int num = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
                this_rate[j] *= omega_kappa_pow[num];
            }
        }
    }
}

void ModelCodon::computeCodonRateMatrix_2KAPPA() {
    int nrates = getNumRateEntries();
    memcpy(rates, empirical_rates, nrates*sizeof(double));

    int i, j;
    double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
    double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};
    double kappa2_pow[] = {1.0, kappa2, kappa2*kappa2, kappa2*kappa2*kappa2};

    for (i = 0; i < num_states; i++) {
        double *this_rate = &rates[i*num_states];
        int *this_rate_attr = &rate_attr[i*num_states];
        if (phylo_tree->aln->isStopCodon(i)) {
            continue;
        }
        for (j = 0; j < num_states; j++) {
            int attr = this_rate_attr[j];
            if (this_rate[j] == 0.0) continue;
            if (attr & CA_SYNONYMOUS) { // synonymous
                int numts = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);            
                int numtv = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
                this_rate[j] *= kappa_pow[numts]*kappa2_pow[numtv];
            } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
                int numts = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);            
                int numtv = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
                this_rate[j] *= omega_kappa_pow[numts]*kappa2_pow[numtv];
            }
        }
    }
}

double ModelCodon::computeEmpiricalOmega() {
    double dn = 0.0, ds = 0.0;
    int i, j;
    if (ignore_state_freq) {
        for (i = 0; i < num_states; i++) {
            if (phylo_tree->aln->isStopCodon(i))
                continue;
            double *this_rate = &rates[i*num_states];
            int *this_rate_attr = &rate_attr[i*num_states];
            for (j = 0; j < num_states; j++)
                if (this_rate_attr[j] & CA_NONSYNONYMOUS)
                    dn += state_freq[i]*this_rate[j];
                else
                    ds += state_freq[i]*this_rate[j];
        }
    } else {
        for (i = 0; i < num_states; i++) {
            if (phylo_tree->aln->isStopCodon(i))
                continue;
            double *this_rate = &rates[i*num_states];
            int *this_rate_attr = &rate_attr[i*num_states];
            for (j = 0; j < num_states; j++)
                if (this_rate_attr[j] & CA_NONSYNONYMOUS)
                    dn += state_freq[i]*state_freq[j]*this_rate[j];
                else
                    ds += state_freq[i]*state_freq[j]*this_rate[j];
        }
    }
    return (dn/ds)*(0.21/0.79);
}
    


bool ModelCodon::getVariables(double *variables) {
	int j;
    bool changed = false;
    if (num_params > 0) {
        j = 1;
        if (!fix_omega) {
            changed |= (omega != variables[j]);
            omega = variables[j++];
        }
        if (!fix_kappa) {
            changed |= (kappa != variables[j]);
            kappa = variables[j++];
        }
        if (!fix_kappa2) {
            changed |= (kappa2 != variables[j]);
            kappa2 = variables[j++];
        }
        ASSERT(j == num_params+1);
    }
	if (freq_type == FREQ_ESTIMATE) {
        // 2015-09-07: relax the sum of state_freq to be 1, this will be done at the end of optimization
		int ndim = getNDim();
		changed |= memcmpcpy(state_freq, variables+(ndim-num_states+2), (num_states-1)*sizeof(double));
//		double sum = 0;
//		for (i = 0; i < num_states-1; i++)
//			sum += state_freq[i];
//		state_freq[num_states-1] = 1.0 - sum;

        // BUG FIX 2015.08.28
//        int nrate = getNDim();
//        if (freq_type == FREQ_ESTIMATE) nrate -= (num_states-1);
//		double sum = 1.0;
////		int i, j;
//		for (i = 1; i < num_states; i++)
//			sum += variables[nrate+i];
//		for (i = 0, j = 1; i < num_states; i++)
//			if (i != highest_freq_state) {
//				state_freq[i] = variables[nrate+j] / sum;
//				j++;
//			}
//		state_freq[highest_freq_state] = 1.0/sum;
	}
    return changed;
}

void ModelCodon::setVariables(double *variables) {
	int j;
	if (num_params > 0) {
        j = 1;
        if (!fix_omega)
            variables[j++] = omega;
        if (!fix_kappa)
            variables[j++] = kappa;
        if (!fix_kappa2)
            variables[j++] = kappa2;
        
		ASSERT(j == num_params+1);
	}
	if (freq_type == FREQ_ESTIMATE) {
        // 2015-09-07: relax the sum of state_freq to be 1, this will be done at the end of optimization
		int ndim = getNDim();
		memcpy(variables+(ndim-num_states+2), state_freq, (num_states-1)*sizeof(double));

        // BUG FIX 2015.08.28
//        int nrate = getNDim();
//        if (freq_type == FREQ_ESTIMATE) nrate -= (num_states-1);
//		int i, j;
//		for (i = 0, j = 1; i < num_states; i++)
//			if (i != highest_freq_state) {
//				variables[nrate+j] = state_freq[i] / state_freq[highest_freq_state];
//				j++;
//			}
	}
}

void ModelCodon::setBounds(double *lower_bound, double *upper_bound, bool *bound_check) {
	int i, ndim = getNDim();

	for (i = 1; i <= ndim; i++) {
		//cout << variables[i] << endl;
		lower_bound[i] = MIN_OMEGA_KAPPA;
		upper_bound[i] = MAX_OMEGA_KAPPA;
		bound_check[i] = false;
	}

	if (freq_type == FREQ_ESTIMATE) {
		for (i = ndim-num_states+2; i <= ndim; i++) {
//            lower_bound[i] = MIN_FREQUENCY/state_freq[highest_freq_state];
//			upper_bound[i] = state_freq[highest_freq_state]/MIN_FREQUENCY;
            lower_bound[i]  = MIN_FREQUENCY;
//            upper_bound[i] = 100.0;
            upper_bound[i] = 1.0;
            bound_check[i] = false;
        }
	}
}

double ModelCodon::optimizeParameters(double gradient_epsilon) {
	int ndim = getNDim();
	
	// return if nothing to be optimized
	if (ndim == 0) return 0.0;
    
	if (verbose_mode >= VB_MAX)
		cout << "Optimizing " << name << " model parameters..." << endl;


	double *variables = new double[ndim+1];
	double *upper_bound = new double[ndim+1];
	double *lower_bound = new double[ndim+1];
	bool *bound_check = new bool[ndim+1];
	double score;

    for (int i = 0; i < num_states; i++)
        if (state_freq[i] > state_freq[highest_freq_state])
            highest_freq_state = i;

	// by BFGS algorithm
	setVariables(variables);
	setBounds(lower_bound, upper_bound, bound_check);
    if (phylo_tree->params->optimize_alg.find("BFGS-B") == string::npos)
        score = -minimizeMultiDimen(variables, ndim, lower_bound, upper_bound, bound_check, max(gradient_epsilon, TOL_RATE));
    else
        score = -L_BFGS_B(ndim, variables+1, lower_bound+1, upper_bound+1, max(gradient_epsilon, TOL_RATE));

	bool changed = getVariables(variables);
    // BQM 2015-09-07: normalize state_freq
	if (freq_type == FREQ_ESTIMATE) { 
        scaleStateFreq(true);
        changed = true;
    }
    if (changed) {
        decomposeRateMatrix();
        phylo_tree->clearAllPartialLH();
        score = phylo_tree->computeLikelihood();
    }
	
	delete [] bound_check;
	delete [] lower_bound;
	delete [] upper_bound;
	delete [] variables;

	return score;
}


void ModelCodon::writeInfo(ostream &out) {
    if (name.find('_') == string::npos)
        out << "Nonsynonymous/synonymous ratio (omega): " << omega << endl;
    else
        out << "Empirical nonsynonymous/synonymous ratio (omega_E): " << computeEmpiricalOmega() << endl;
    out << "Transition/transversion ratio (kappa): " << kappa << endl;
    if (codon_kappa_style == CK_TWO_KAPPA) 
        out << "Transition/transversion ratio 2 (kappa2): " << kappa2 << endl;
}

